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^ ! Abstract 

^ \ In this manuscript, we analytically and numerically study statistical properties of an het- 
\ eroskedastic process based on the celebrated ARCH generator of random variables whose 
variance is defined by a memory of g m -exponencial, form (e® m=1 = e x ). Specifically, we 
inspect the self-correlation function of squared random variables as well as the kurtosis. 
In addition, by numerical procedures, we infer the stationary probability density func- 



tion of both of the heteroskedastic random variables and the variance, the multiscaling 

properties, the first-passage times distribution, and the dependence degree. Finally, we 

introduce an asymmetric variance version of the model that enables us to reproduce the 

c/j ■ so-called leverage effect in financial markets. 
O 

1 Introduction 

Many of the so-called complex systems are characterised by having time series with a 
peculiar feature: although the quantity under measurement presents an autocorrelation 
function at noise level for all time lags, when the autocorrelation of the magnitudes is 
appraised, a slow and asymptotic power-law decay is found. This occurs, e.g., with (log) 
price fluctuations of several securities traded in financial markets p], temperature fluctua- 
tions [2], neuromuscular activation signals [3] or even fluctuations in presidential approval 



> 



O ■ ratings jl] amongst many others [S]. Moreover, most of these time series are also char- 
acterised by probability density functions with asymptotic power-law decay and a profile 
suggestive of intermittency that is identified by regions of quasi-laminarity interrupted 
by spikes. In this perspective, this type of time series might be seen as a succession of 
measurements with a time-dependent standard deviation. Mathematically, this type of 
stochastic process is defined as heteroskedastic, in opposition to the class of processes with 
constant standard deviation that is defined as homoskedastic. With the primary goal of 
reproducing and forecasting inflation time series, it was introduced in 1982 the autore- 
gressive conditional heteroskedasticity process (ARCH) [6]. The ARCH process rapidly 
has come to be a landmark in econometrics giving raise to several generalisations and 
widespread applications not only in Economics and Finance but in several other fields as 
well. 

In the sequel of this article, we introduce further insight into a variation of the ARCH 
process studied in Ref. [7] which is able to reproduce the properties we have referred to here 
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above. Our considerations are made both on analytical and numerical grounds. Although 
the primary goal of this manuscript is an extensive description of the model following 
the lines of Ref. [7], some assessment of its capability of reproducing the same features 
of SP500 daily log fluctuations spanning the 3 rd January 1950 to the 28 th February 2007 
is made. In this context, we also introduce a slight modification on the model which 
turns it able to reproduce the leverage effect when the model is applied to surrogate price 
fluctuations time series. The manuscript is organised as follows: after introducing the 
ARCH processes and present some general properties, we make known in Sec. [2] some 
analytical calculations on the autocorrelation function of the model herein analysed and 
the correlation between variables and squared standard deviation for the extension as well 
as the expressions for the kurtosis. In Sec. El we introduce results from the numerical 
analysis about the probability density functions of the stochastic variable, z t , and its 
squared instantaneous standard deviation, of; The dependence degree between z\ and 
2f +r ; The distribution of first passage times of ; and the multiscaling properties. In 
Sec. HI we establish an asymmetric variation of the model which allows the reproduction 
of the so-called leverage effect. Final considerations are addressed to Sec. El 



2 The symmetric variance model 

We start defining an autoregressive conditional heteroskedastic (ARCH) time series as a 
discrete stochastic process, z t , 

Zt = a t u t , (1) 

where ujt is an independent and identically distributed random variable with mean equal 
to zero and second-order moment equal to one, i.e., (u t ) = and (of) = 1. Usually, uj is 
associated with a Gaussian distribution, but other distributions of uj have been presented 
to mainly describe price fluctuations [8J. In his benchmark article of Ref. [Bj, Engle 
suggested a dynamics for of establishing it as a linear function of past squared values of 

s 

a* = a + J2 b izLn M*>0). (2) 
i=i 

In financial practise, namely price fluctuation modelling, the case s — 1 (b\ = b) is, by far, 
the most studied and applied of all ARC HASke processes. It can be easily verified, even 
for all s, that, although (z t z t >) ~ 8 tt ', correlation (\z t \ \zp\) is not proportional to 5 tt '- 
As a matter of fact, for s — 1, it has been proved that, (zf z#) decays as an exponential 
law with a characteristic time r = |ln&|~~ , which does not reproduce empirical evidences. 
In addition, the introduction of a large value for parameter s bears implementation prob- 
lems [10J. Expressly, large values of s soar the complexity of finding the appropriate set of 
parameters {bi} for the problem under study as it corresponds to the evaluation of a large 
number of fitting parameters. Aiming to solve the imperfectness of the original ARCH 
process, the GARCH (s,r) process was introduced [H], with Eq. (j2J) being replaced by, 

s r 
i=l i=l 
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Nonetheless, even this process, presents a exponential decay for (z* z 2 ,), with r = |ln (b + c) 
for GARCH (1,1), though condition, b + c < 1, guarantees that GARCH (1,1) corre- 
sponds exactly to an infinite-order ARCH process |12j . 

Despite the fluctuation of the instantaneous volatility, the ARCH(1) process is actu- 
ally stationary with a stationary variance, given by, 

<- 2 > = - 2 = I ^, (6>1), (4) 

((...) represents averages over samples and ^T. averages over time). Moreover, it presents 
a stationary probability density function (PDF), P(z), with larger kurtosis than the 
kurtosis of distribution P(oS). The fourth-order moment is, 



This kurtosis excess is precisely the outcome of time-dependence of a t . Correspondingly, 
when 6 = 0, the process reduces to generating a signal with the same PDF of uo, but with 
a standard variation y/a. 

In the remaining of this article we consider a ARCH(1) process where an effective 
immediate past return, z t -i, is assumed in the evaluation of of [7j. Explicitly, Eq. (T2]) is 
replaced by 

a 2 t =a + bz 2 _ l , (a,b>0), (5) 
where the effective past return is calculated according to 

t-i 

X)£(i-t + l)z?, (6) 



t-i 

= 

i=to 



where 



T 



(f < 0, T > 0, q m < 2) (7) 



with 

exp g [x] = e* = [1 + (1 - q) x]f~* , (q G 7Z) , (8) 

Z qm (t') ee Y^H=-t' ex Pg m [t] (N+ = max {0' known in the literature as q-exponential 

This prososal can be enclosed in the fractionally integrated class of heteroskedastic pro- 
cess (F I ARCH). Although it is similar to other proposals [Hj, it has a simpler structure 
which permits some analytical considerations without introducing any underperformance 
when used for mimicry purposes. For q = — oo, we obtain the regular ARCH (1) and for 
q — 1, we have /C (f) with an exponential form since exp x [x] = e x [T5]. Although it has a 
non-normalisable kernel, let us refer that the value q m = oo corresponds to the situation 
that all past values of Zt have the same weight, 



2 This condition is known in the literature as Tsallis cut off at x = ± (1 — q) . 
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In this case, memory effects are the strongest possible, i.e., every single element of that 
past has the same degree of influence on af making it constant after a few steps. Because 
of this, in this case, P(z) is the same as noise uj, as shown in [7J. Similar heuristic 
arguments are the base for the Gaussian nature of the "Elephant random walk" [16] . 

Assuming stationarity in the process some calculations can be mad^f). Namely, it is 
provable that the average value of a 2 yields, 

(a 2 ) = a* = ^ ( 6> i), (10) 

and the covariance, (z t z t >) corresponds to 

(z t z t >) = (a t a v ) (u t LO t >) , (11) 

which, due to the uncorrelated nature of u, gives (z t z t >) = for every t ^ t' and (zf) = 
(a 2 ). In addition, we can verify that all odd moments of z t are equal to zero. Concerning 
the fourth-order moment, (zf), we have 

(zf) = (afafu J 2 t u J f) = ([af] 2 )(u J f), (12) 

which by expansion yields, 

(zf) = a 2 (oof) + 2ab^£K(i-t + l) (zf) (ujf) + 

i=to 

b 2 Z[K:(i-t + l)] 2 (zi)(ut)+ (13) 

i=to 

2b 2 E E )C(i-t + l))C(j-t + l)(z 2 z 2 )(ujf). 

i=to j=i+l 

If Zi and Zj are assumed as strictly independent, then {zf z 2 ) = (zf) (z 2 ) = (zf) 2 . 
Assuming stationarity we have, 

2 /, ,4\ 



[zf)= a 2 (uf) + 2ab(zf)(uf) + 

b 2 (zf) (ujf) E [K{i-t+ l)f 

l=to 



2b 2 (zf) 2 (uf)Z E K(i-t + l)K.(j-t + l) 



(14) 



or 



r 4 N a 2 + 2ab(zf) + 2b 2 (zf) 2 Q, , 4N 

h) i~ 1 - b 2 (uof) Q 2 (15) 



3 This has been numerically analysed by computing the self-correlation function for different waiting 
times 
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with Q x = £ E K(i-t + l)K.{j-t + 1) and Q 2 = £ [K(i-t+ l)] 2 . On the 

«=*o j=i+l i=to 

other hand, we have the other limiting case, (zf 2J) = ). Equation (Tl3l) is then written 

as 

^^-l-6M^)(2Qx + Q 2 ) H> - ( 6) 

The labelling as upper bound for (zf)j and lower bound for (zf) c comes as follows; The 
introduction of non-Gaussianity in heteroskedastic processes comes from the fluctuations 
in the variance (or in zf), when the variables are strongly attached between them, there 
is a small level of fluctuation in of and eventually it becomes constant. With <j t being 
a constant, or approximately that, there is not introduction of a significant level of non- 
Gaussianity measured from (zf), hence (zf)j > (zf) c (o^O). 

For an accurate description of (zf), which lies between the two limiting expressions, 
we must compute correlations (zf zf). That is obtained averaging, 

t-i 

zf zf = a 2 uf uof + a b £ (* — t + 1) zf ujf uf + 

i=to 

ab t! ZJC(j-t' + l) z]ufuf+ (17) 

3=to 

b 2 E E K (i ~ t + 1) K (j - t> + 1) zf z) uf uf . 

i=to j=to 

Defining r = t — t', the last term of rhs of Eq. (1171) . hereon labelled as C, is responsible 
for the dependence of (zf zf) with r. It can be written as 

C = b 2 uf ajj § *~E 1 JC(i-t + l)K{j -t + l + r) zf z] 

*=*o j=to 

f E JC(i-t + l)1C(i-t + l + T) zf+ (18) 

i=to 

E _1 E K,(i-t + l + r)lC{j-t + l)z 2 z 2 . 

i=to j=t—T+l 

We shall now consider a continuous approximation where the summations are changed by 
integrals, 

. . . — > / ... dx, 

and t ^> 1, so that the following relations are obtained in the limit t — > oo. Computing 
(C), the term in zf has the coefficient, 

Ci (r) ~ ex P(?m [x - t) exp 9m [x - t + r] dx 

(19) 

2-q-n 



t^OD 



[r(q m -l)]l-gm p 

2-q m 2 > l 



— ■ 2 + 

9m— 1' 9m— 1' ?m-l' (l-<?m)r 



where F 2y i [a, b; c; x] is the generalised hypergeometric function. Regarding the term in 
z 2 z 2 , the first approximation is obtained considering (zf z|) = (zf) = (z 2 ) 2 . Its 
coefficient is then given by 

-t-T 



c 2 (r) ~ / t _ T J T exp 9m [y - t] ex P(/m [a; - t + r] rfa; 

exp 9c [-A r] , 



(20) 



t— >oo 



with 



g c = (21) 

and A = q~ x . A simple inspection shows that C 2 (r) decays much slower than C\ (t), 
hence the asymptotic form of {z 2 z 2 ) is dominated by C 2 (r) as it is illustrated in the inset 
of Fig. [U In Fig. (TJ we bring face to face Eq. (1201 and the autocorrelation function of 
z 2 from numerical simulation using the parameters applied to reproduce SP500 returns 
previously determined, namely {a, b = b SP ,q m = qsp} = {0.5,0.99635, 1.6875}. 

From all these equations we are able to conjecture expressions which relate parameters 
{a, b, q m } with the form of the distribution for the case where (zf) is finite. In this way, 
we can use the ansatz that the distribution of this dynamical model is associated with a 
g-Gaussian (or Student-t) distribution @, 

p(z)=Aexp q [-Bz 2 ], (g<3), (22) 
with B = [a 2 (3 — q)] 1 , where, 

a\ = J z 2 [ P {z)} q dz/ J \p{z)fdz, 

is the g-generalised second order moment [13], and A is the normalisation constant. This 
assumption is based on the same type of arguments used in [27] and whose accuracy we 
verify later on (see Sec. [3]). For q < 5/3, o 2 q relates to the usual variance according to 
a 2 (3-g) = a 2 (5-3g). 

From Eqs. (jI5), 071), (fl9"|h and (120) we can write, 



(zf) = a 2 (uf) + 2a b (z 2 ) (uf) + b 2 (zf) (uf) ^ 



+b 2 (u?) a 2 (to 2 ) 2 + 2ab (z 2 ) (to 2 ) 2 + 2/C 



(23) 



where K represents terms like, 

/ / ex P<? m l x ~ A ex Pg m [x + r-t] [Ci (r) + C 2 (r)] dr dx, 
Jo Jo 



(24) 



4 A g-Gaussian, with q > 1, corresponds to a Student-t with m degrees of freedom with q = where 
m is taken as a real positive number. 
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which corresponds to a quite complex integration over x of hypergeometric function special 



functions like F' 



2.1 



c represent genera' 



a, b; c; x 
values. 



and the Appell hypergeometric function [T7j where a, b, and 



Taking into attention that for a g-Gaussian with q < ~ its fourth moment is, 

2 3 g — 5 



<^> = 3(a 2 )- 



5g-7' 



(25) 



we can obtain approximate relations between the parameters of the model and the param- 
eters of the distribution. This is achieved when we equalise Eqs. ( |23|) and ( |25i) . remember- 
ing the expression of the variance, Eq. ( flOi) . and the form of the autocorrelation function 



of z 2 . This procedure is obviously important in parameter estimation. Therefore, from 
the decay of the z 2 autocorrelation function we can determine the value of q m , and a and 
b from the equalisation we have just referred to together with Eq. (TiT)]) . 



3 Numerical considerations 



3.1 Stationary probability density functions 

We firstly recover the results previously presented for the adjustment of P (z) with q- 
Gaussians. As mentioned in [7] the adjustment is striking (diagrams of q as a function of 
b and q m are presented in Figure 1 of that reference). Using the method of x 2 minimisation, 
we have obtained for the same cases previously studied^ average values of x 2 — 1.1 x 1CT 6 
(per degree of freedom) and R 2 = 0.99990. In Fig. [2j we present an example for which it 
is possible to assent the accuracy of the fitting not only in the tails, but in the peak of 
the PDF as well. We have performed further analysis using the cumulative distribution 
function (CDF) and the Kolmogorov-Smirnov Distance, D^Si 

D KS = max\H{z)-H {z)\, (26) 

where H (z) is the empirical CDF obtained from numerical evaluation of the model and 
H (z) is the testing probability density function, 

H (z)= I A exp q [-Bx 2 ] dx. (27) 

J — oo 

The average Dks value obtained for the same cases plotted in Figure 1 of the prior work 
is equal to 4.25 x 10 -3 . Such values allow us to rely on the null hypothesis [T9] . 

p ( z ) = p ( z ) = A ex P(? [-B z 2 ] . (28) 

Based on the acceptance of the null hypothesis (|28p we are able to introduce some 
insight into the distribution of cr 2 , p a (a 2 ). Firstly, we carry out the following change of 
variables, 

z t = In z 2 

a t = In a 2 , (29) 
Ut = In oj 2 

5 Comparing with prior studies we have increased the runs by a factor of 10. 
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so that Eq. (TjQ) turns into, 

z t = a t + u t . (30) 
In probability space, regarding that a t and ui t are independent, we have, 



p{z) 



Pa ip) Pu, (z - da. 



(31) 



We can now apply the convolution theorem. Being p (z t ), the probability of z t , according 
to such a theorem, 

p(z t ) =T- 1 [P»{k)P^{k)] , (32) 

where 



P x (k) 



2 7T 



/ (x) exp [ikx] dx = T [f (x)) 



(33) 



and T 1 [f x (k)] = -4= J f x (k) exp [—i kx] dx = f (x) is the inverse Fourier Transform. 
Since we respectively know and postulate the form of p (u) and p (z), we can write down, 



p{u) 



/2tt 



exp 



2 2 



(34) 



p(z) =A exp q [-Be s ] exp [§] , 
(B = B (a = 1)), yielding the respective Fourier Transforms [20J, 



+ 



(35) 



r[i + z/:]F 2 ,i -ri,| + iA;,| + ^,i3-i3g 



where Q = jz^ + ^ + ik, B [. . .] is the Beta function, and F 2i i [. . .] is the regularised hyper- 
geometric function [17J. Applying Eq. (1351) in Eq. ( 1321) we can compute the distribution 
of <7 (easily related to p a (cr)), 



P* (a) = 



A(fc) 
Pa(k) 



From a laborious and tricky calculation, using properties of P2,i 
properties), it can be verified P„ (a) cooresponds to, 



T{9-ik] 

(2f3) lk T[c} 



(36) 

(see [18J and related 
(37) 



Therefore, a 2 follows an inverse Gamma distribution, 

1 



I 2\- 1 ~ c 

[a J exp 



(2^) c r[c] 



26a 2 



(38) 



where c = and 9 = ^^} 3q ^ ■ This result attests the validity of the superstatistical 
approach to the problem of heteroskedasticity. It is worth mentioning that superstatistics 
[2T] represents the long-term statistics in systems with fluctuations in some characteris- 
tic intensive parameter of the problem like the dissipation rate in Lagragian turbulent 
fluids [22] or the standard deviation like in the subject matter of heteroskedasticity. For 
the values of q$p and b$p, with have obtained random variables z t associated with a 
(q = 1.465)-Gaussian. This yields csp = 1.648 . . . and 9sp = 0.770 . . ., which have been 
applied in Fig. [3] to fit p a (er 2 ) obtained by numerical procedures. In that plot, it is visible 
that numerical and analytical curves are in proximity. 



3.2 Dependence degree 

The degree that the elements of a time series are tied-in is not completely expressed 
by the correlation function in the majority of the cases. In fact, regarding its intimate 
relation with the covariance, the correlation function is only a measure of linear depen- 
dences. Aiming to assess non-linear dependences, information measures have been widely 
applied [23] • In our case, we use a non-extensive generalisation of Kullback-Leibler infor- 
mation measure 24, 25M, 



h = ~^p{z 2 t) z 2 t+T ) 



p{ z h z t+ T ) 



l-q 



(39) 



where p' (z 2 ,z 2 +r ) = p(z 2 )p(z 2 +T ) = [p(z 2 )] 2 (assuming stationarity), which has been 
able to provide a set of interesting results with respect to dependence problems [2B] ■ The 
quantification of the dependence degree is made through a value, q op , which corresponds 
to the inflexion point of the normalised version of I q , 

R q = (40) 

where 7™ ax is the value of I q when variables z (t) and z (t + r) present a biunivocal de- 
pendence (see full expression in Ref. [25]). 

For infinite signals it can be shown that, when the variables are completely independent 
q°p = oo, whereas q op = when variables are one-to-one dependent. For finite systems, 
there is a noise level, which is achieved after a finite time lag r. Typical curves of R q 
are depicted in Fig. [5] for q m = qsp and b = bsp- 

In what follows, we present results obtained from numerical evaluation of R q for differ- 
ent values of r. As expected, dependence relies on the balance between the extension of 
memory, which is given by q m and the weight of effective past value, Zt-i, on o t . Firstly, let 
us compare cases {q m = q$p, b = b SP } and {q m = 1.25, b = b SP }, as an example of what 
happens when we fix b as a constant (see Fig. H]). Dependence is obviously long-lasting 
in the former case, in the sense that it takes longer to attain g° p , but for small values of 



3 In the limit q — > 1 the Kullback-Leibler mutual information definition is recovered. 
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r, the latter has presented higher levels of dependence. This has to do that /C (£') is nor- 
malised and that implies the intersection of the curves for different values of q m at some 
value of t' . Alternatively, when q m decreases, the recent values of z t have more influence 
on Zt-i than past values. When the value of q m is kept constant, we have verified that 
smaller values of 6 lead to a faster approach to noise value q op . In a previous work on 
GARCH (1,1) [27J, we have verified that variables approximately associated with the 
same distribution present the same level of dependence independently of the pair (6, c) cho- 
sen. In this case, recurring to cases {q m = 1.375, 6 = 0.75} and {q m = 1.625, b = 0.875}, 
we have noticed that the curves present very close values for small lags, but they fall apart 
for r > 10, revealing a more intricate relation between q, q m , and b than in GARCH (1, 1). 
Additionally, comparing dependence and correlation, we have verified that the decay is 
faster for the latter. Specifically, taking into account noise values of q op and C T (z^), it is 
verifiable that q op takes longer to achieve q op than C T (zf ) takes to reach C n (z%). 



3.3 First-Passage Times 

First-passage studies in stochastic processes are of considerable interest. Not only from the 
scientific point of view [2H] (it is useful in the approximate calculation of the lifetimes of 
the problems/systems) as well as from a practical perspective, since they can be applied 
to quantify the extent of reliability of forecasting procedures, e.g., in meteorology or 
finance [U [29j [30]. In what it is next to come, we have analysed the probability of 
zf G Si = [a,b) and z^ +T e Si. We have divided zf domain into five different intervals. 
Explicitly: 



• 




A < i; 


• 


S 2 : 


1< z 2 t < 2; 


• 


S 3 : 


2<z 2 t < 5; 


• 


S 4 : 


5 < z\ < 10 


• 


S 5 : 


zf > 10. 



Analysing the probability density functions we have verified that the simplest expres- 
sion which enables a numerical description of first-passage inverse cumulative probability 
distribution, T> (t), is a linear composition of a asymptotic power- law (or a //-exponential) 
with a stretched exponential, 



V(t) 



l + (v 



t 



1 

1-1/ 



+ (1 - e) exp 



t 

^2 



(41) 



Curves of some analysed cases are presented in Fig. [U] and fitting parameters in Tab. HJ The 
cases we present are: I-{q m = qsp, b = 0.5}, II-{q m = 1.25, b = b S p}, III-{q m = 1.25, 6 = 0. 
IV-{q m = 1.5,6 = 0.875}. 

From the Fig. [6] we have verified that, excepting region Si, all curves of T> (t) exhibit 
a decay closely exponential [y — 1). For region Si, as we increase the non-Gaussianity of 
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p(z t ), we have observed that both of the values of e and v approach one. Comparing the 
remaining curves, we have verified that, for every region S2 — 5*5, the set of parameters 
which provides higher degree of non-Gaussianity has the larger characteristic times X. 
Keeping the memory parameter q m constant, we have observed that the higher b, the 
higher X is. An inverse dependence is found when we have fixed b and let q m vary. 
In other words, smaller values of q m (which enhance broader distributions) have larger 
values of X. We have also verified that the first-passage times are not Poisson distributed 
as it is straightforwardly verified in S% plots. Looking at the values of squared daily 
index fluctuations of SP500, we have verified that they present an asymptotic power-law 
behaviour for T> (t) (e = 1, v > 1) (fitting parameters shown in Tab. [2]). Comparing the 
results of the model with empirical results from SP500 time series, we have observed that 
the model provides an overall reasonable description of first-passage times with curves 
almost superposing for S 2 — S4 regions. For Si and S5 regions curves present similar 
exponents, but different values of X. It is worth remembering that we can improve the 
results by considering some characteristic time in K, (t) equation. Taking into account 
the v exponents obtained for the adjustment of 5P500 first-passage times, we verify that 
larger and smaller v- exponents are quite different. Such a strong gap invalidates, at 
a daily scale, the collapse (existence of a single exponent) of V (t) curves proposed for 
high-frequency data [29] . 



3.4 Multiscaling properties 

Multiscaling has been the focus of several studies in the field of complexity [31], par- 
ticularly regarding applications to finance [32j[33j[34j. If in many works multiscaling 
(multifractal) properties of price fluctuations are presented, other studies have put those 
multiscaling properties into questioning (35]. In this section, we analyse mainly multi- 



fractal properties of z t and z+. To this aim, we define a generic variable 



,(2) 



,(2) 

't+i-1' 



(42) 



i=l 



From it, we compute, 



n h (T) = ( Zf ] (T) - (Zf> (T) 



r(2) 



If there is multiscaling, then the following scaling property is observed, 

n h (T) oc T^ h) . 



(43) 



(44) 



The computation of Qh (T) has been made using the well-known MF-DFA procedure [36] . 

For the case of z t , the multiscaling can be easily and analytically explained. The 
multiscaling properties of a time series can emerge twyfold: from memory and from non- 
Gaussianity. Since, by definition, the heteroskedastic process we present is uncorrelated, 
then the only contribution to multiscaling comes from the non-Gaussian character of 
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probability density functions. In this way, time series {z t } is not a multifractal, but a 
bifractal instead [SB]. Therefore, the r\ (h) curve is defined as follows, 



SP500 time series and the surrogate data have enabled us to verify that the model is 
adequate to reproduce the scaling properties of Zf which are basically linear according to 
r) (h) = h. For a constant b value, we have observed that higher correlations, introduced 
by increasing the value of q m turn multiscaling properties weaker. In other words, rj z 2 (h) 
approaches the straight line rj z 2 (h) = h. Similar qualitative results have been presented 
for traded value where long-lasting correlations dominate specially for highly liquid equi- 
ties [34] . Freezing memory parameter, q m , we have verified similar results, i.e., increasing 
the value of b, we increase the tails in zf forcing the multiscaling curve to divert (even 
more) from the straight line t] = h. The same effect is obtained when memory is short- 
ened by reducing the value of q m . By this we mean that, when memory decays faster we 
have a detour from the straight line 7] = h and an approximation to a bifractal behaviour 
because of the asymptotic power-law behaviour. This reflects the fact that the dynamical 
and statistical properties of our system strongly depend on the "force relation" between b 
and q m . Some results from which these observations can be confirmed are shown in Fig. [7J 

4 Asymmetric variance model 

In several systems it has been verified that the correlation function between the observable 
and its instantaneous variance exhibits an anticorrelation dependence. For example, this 
occurs in the case of financial markets, when the correlation between past price fluctua- 
tions and present volatilities is measured. The shape of the curve copes with the so called 
leverage effect [371 EH]. This feature is intimately related to the risk aversion phenomenon, 
i.e., falls in price turn the market more volatile than price rises. In order to reproduce 
this characteristic we introduce a small change in Eq. ©, specifically, 



where c > 0. It is worth stressing that this modification does not introduce any skewness 
in the distribution P (z) which is still symmetric. It only acts on how positive and negative 
values of z t , with the same magnitude, influence z^ by different amounts. 

Using Eqs. (JSJ) and fj4T)]) in Eq. ([T]) and expanding it up to first order we have 





(46) 




(47) 



k i=to ) 

Computing z t zJ +T the only terms which do not vanish after averaging are 



t-l+T 



7i = -2y/a~bc £ (i - (t - 1 + r)) z\ uj t col 




i=t 
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and 

T 2 = -26c 3 * if E)C(i-(t-l + r)) /C (j - t + 1) x 

i=to j=t r^g\ 
zfzjujt UJ 2 t+T . 

Performing averages and considering stationarity we have in the continuous limit, 



(7!)= --6y^bc£ +T JC(x-(t + r)) 
(z 2 ) (a t ) (uf) 5{x-t) dx, 

which gives, 



x 



(50) 



(7i> ~ - exp gm [-r] . (51) 
T 2 = -26c 3 | j'l K (x - (t + r)) K,{x-t) z 7 x uo t uj 2 +T dx+ 

2 £ £ £ (x - (t + t)) /C (y - ^ 3 4 tot u 2 t+T dx dy+ (52) 

It +T lo JC(x-(t + T))JC(y-t) z\z^ u t u 2 +T dxdy}. 
Averaging, only the first integral has a non-null contribution yielding, 

(T 2 ) = -2 6 c 3 (2 - q m fex Vqm [-r] {z Q t ) (a t ) {uj 2 ) . (53) 
It is not hard to show thatQ 

(44^)= 5 6 (x - t) (a t ) (z 2 ) (zf) + 

2A6(y-t)5(x-y) (a t ) (a 2 ) (z 2 ) 2 + 

72 5(x-t)[5(x-y)] 2 (a t )(a 2 )(z 2 ) + (54) 

12 5{x-y) (a t ) (z 2 x z 3 y u t ) + 

12 5(y-t)5(x-y) (a t ) 2 (z 2 x z 2 ) . 

Inserting Eq. (I54p in Eq. (1521) the last two integrals give zero. Therefore, in the first 
approximation, the leverage function, 

goes as, 

£~ (%) + (%) ~-exj>.[-T]. (56) 



7 In order to keep the formulae as simple as possible we use the following expressions the discrete 
notation z t . Formally, it should be read as z (t) since we are dealing with a continuous approach. 
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As it can be seen from Fig. [H] the approximation provides a satisfactory approximation 
of the numerical results. A precise description can obviously be obtained by considering 
higher-order (slowly decaying) terms which are obtained through a quite tedious compu- 
tation that follows exactly the same lines we have just introduced. 

The result above is in apparent contradiction with previous work in which an expo- 
nential dependence with r is defended in lieu of an asymptotic power-law dependence. 
Nevertheless, in Fig. we show the leverage function computed from the 5P500 and 
numerical adjustments with function, 



with q m = q S p and q m = 1. Computing the adjustment error, \ 2 an d R 2 , we verify that 
both approaches present similar values for the numerical adjustment with q = qsp being 
scanty better. For q = q$p, we have obtained x 2 = 4 x 1CT 5 and R 2 = 0.47 and for the 
exponential adjustment x 2 = 5 x 10~ 5 and R 2 = 0.44. From these results, we can affirm 
that our proposal is, at least, as good as the exponential decay scenario firstly introduced 
in Ref. [3B]- It is also worth noting that, although this variation is asymmetric concerning 
the effects of the sign of variables z t on the evaluation of a 2 , the model is ineffective in 
reproducing the skewness of the distribution of price fluctuations. This owes to the fact 
that Ut used up to now is symmetric thus, also annulling the asymmetry introduced by 
the Eq. (146p . It is worth stressing that moments (z^), for n even, in this section are 
obviously different from the values presented in preceding sections. 

5 Final remarks 

In this manuscript we have introduced further insight into a heteroskedastic process en- 
closed in the class of fractionally integrated ARCH processes. This process is charac- 
terised by a memory of past values of the squared variable which decays according to an 
g m -exponential. Despite of the fact that we were unable to provide an analytical proof, 
prevailing statistical testing has shown that g-Gaussian distributions properly describe 
the probability density function of the generated stochastic variable. Based on this fact, 
we have determined a form of the instantaneous variance, a 2 , probability density function 
which has yielded a inverse Gamma distribution like it happens in superstatistical models 
giving g-Gaussians as long-term distributions. Moreover, we have computed the first term 
of the correlation function, C T (z 2 ), which corresponds to a g c -exponential. An analytical 
relation between q m and q c is presented. From these results, we are able to state that this 
dynamical system can actually be described within non-extensive statistical mechanics 
(NESM) framework by a triplet of q values |39j. As a matter of fact, this process presents 
all the elements to be characterised, in z 2 variable, as a NESM process. Explicitly, besides 
presenting asymptotic power-law distributions which maximise non-additive entropy S q 
(as z t ), it has a slow decaying (g c -exponential) auto-correlation function, and it exhibits 
multiscaling properties. Such properties have been advocated as primary features of sys- 
tems that should be studied within NESM framework (for related literature see [3D]) for 
a long time. Furthermore, we have verified that, for a sufficiently high level of mem- 
ory the model presents a non-exponential distribution of first-passage times and strong 




(57) 
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levels of dependence measured from a generalisation of Kullback-Liebler mutual infor- 
mation. Though they have been applied in several other areas, in view of the fact that 
hetoskedastic processes have been introduced in a financial context, we have tested the 
model against daily index fluctuations of SP500. The results firstly presented, together 
with the results of this manuscript, show that this model, despite of its simplicity (it only 
has two parameters), is able to reproduce the most relevant and important properties, 
namely the probability density functions, the Hurst exponent, the autocorrelation func- 
tions, multiscaling, and first-passage time distribution (in a less good extend compared 
to previous). It is our belief that the same occurs with other time series presenting sim- 
ilar characteristics. Moreover, we have studied how the extension of the memory (tuned 
by q m ) and its weight, or alternatively, the weight of fluctuations in of (adjusted by b) 
have on the quantities. Still in the context of financial time series, we have introduced 
a slight modification which allows the reproduction of the leverage effect. Under these 
circumstances, we propose that the leverage is not described for an exponential function, 
but for a g m -exponential function instead. When statistically tested, this proposal has 
emerged as good as the exponential description. It is well-known that many distributions 
obtained from complex systems present skweness which this model has not been able to 
capture because of the symmetrical nature of noise u> distribution. The use of other types 
of noise u, jointly with the modification presented in Sec. H] might give rise to even more 
precise modelling. 
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Figure 1: The black line presents the numerical evaluation of the z\ autocorrela- 



tion function vs the lag r and the red line C2 (r) from Eq. (1201) for parameters 
{0.5,6 = bsp,q m = Qsp}. The inset depicts the way C\ (r) (dotted line) and C2 (r) (full 
line) decay. As it can be seen, C\ (r) decays much faster than C 2 (t) which is the major 
responsible for C T (zf) behaviour for large r. 
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Figure 2: Upper panel: The symbols represent the probability density function P (z) vs 
z for parameters q m = 1.5 and b = 0.875, and the line the best fit for a g-Gaussian with 
q = 1.385 [on log-linear scale] and unitary standard deviation (x 2 = 1.39 x 10~ 6 and 
R 2 = 0.99987). The inset is the same, but on linear-linear scale permitting the appraisal 
of the fitting in the central region. Lower panel: The symbols are for the empirical 
cumulative distribution function H (z) vs z for the same parameters and the line the 
CDF of a g-Gaussian with q = 1.385 and unitary standard deviation (D^s — 0.00457). 
The inset is the same plot, but on a log-linear scale. 
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Figure 3: The black symbols represent p a (a 2 ) vs a 2 [on log-log scale] obtained by numeri- 
cal evaluation of the process with q m = Qsp an d b = bsp yielding a g-Gaussian distribution 
with q = 1.465 (x 2 = 1.6 x 10 -6 ) shown in the inset [on log-linear scale]. The red line in 
the main plot is the representation of the inverted Gamma distribution with c = 1.648 . . . 
and 6 = 0.770 . . . (x 2 = 6.1 X lO" 4 ). 
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Figure 5: Values of q op vs time lag r for the pairs {q m = qsp,b = bsp} as indicated in the 
figure. The dotted lines are merely presented as a guide to the eye. 
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Figure 6: Inverse cumulative distribution T> (t) of the first-passage times vst for the sets 
of parameters mentioned in the text. Upper panels: region S\ (left) and S2 (right); 
Middle panels: region 53 (left) and S4 (right); Lower left panel: region S 5 (left). For the 
regions S 2 — 5*5, the exponential decay is evident. Lower right panel: Inverse cumulative 
distribution V (t) of the first-passage times vs t for squared daily price fluctuations of 
SP500 (symbols) with the best fit represented by the dotted lines. The full lines are 
obtained from the model with parameters q m = qsp and b = b SP . The non-exponential 
behaviour is clear in this latter case. The curves of the S3 — S5 regions are shifted by a 
factor of 2, 4, and 16, respectively. 
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Figure 7: Multiscaling exponents 77 ws moments h of the values presented in the plot. The 
dashed line is 77 = h. 
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Figure 8: The symbols represent leverage £ vs time lag r of a processes using Eq. (HS1) 
with q m = 1.65, b = 0.95 and c = 0.1 obtained from a time series of 10 6 elements. The 
red line is the best fit using Eq. ([561) with £ (r = 0) = -0.39 (x 2 = 3 x 10" 4 ). 
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Figure 9: The symbols represent the absolute value of leverage £ vs time lag r for the 
daily index fluctuations of SP500 spanning the period mentioned above. The red line 
represents the fitting for a g-exponential function with L = 0.07, q = qsp and T = 9, 
and the green line a exponential fitting with Lq = 0.06 and T = 20 (fitting error values 
in the text). 
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